常见的群落相似性或距离测度的计算

常见群落相似性或距离测度的计算

前篇已经初步讲述了关于群落多样性分析中的Beta多样性基本概念，并提到常见的Beta多样性分析方法一般建立在群落相似性/或距离测度的基础上，以及初步介绍了有关群落相似性/或距离测度的基础。

相似性或距离的衡量标准有很多种，Legendre和Legendre（1998）列出大约30种方法，并对生态相似性作了更详细的介绍，有兴趣可自行参阅Legendre和Legendre（1998）“Numerical Ecology”第七章“Ecological resemblance”的内容。

接下来的篇幅，白鱼小编选择了其中几种常见的相似性或距离测度为大家作简介。前半部分介绍公式，后半部分展示在R语言中的计算方法。

常见的相似性或距离简介

定性（二元）对称相似性指数

这种类型的相似性指数不适用于群落数据分析，忽略。

定量对称相似性指数

同样地，这种类型的相似性指数也不适用于群落数据分析，忽略。

定性（二元）非对称相似性指数

Jaccard相似性/相异度

Jaccard相似性指数（Jaccard similarity index）将两个样方共享的物种数量（a）除以两个样方中出现的所有物种的总和（a + b + c，其中b和c是仅在第一个和第二个样方中出现的物种数量）。计算公式如下：

转换为Jaccard相异度（Jaccard dissimilarity）：

注：如前篇中所提到的，相似性指数（S）和距离测度（D）常见3种转化公式：

D = 1-S，或D = (1-S)^1/2，或D = (1-S²)^1/2；

其中“D=1-S”更为常见；对于“D = (1-S)^1/2，或D = (1-S²)^1/2”类型，这种转化的目的使某些距离指数具有欧式几何特征，欧式属性的距离在某些分析中将会非常有用。这个将来我会在一些更具体的方法中提到。

本篇中所展示的相互转换，默认以“D = 1-S”为例，下同。

Sørensen相似性/相异度

与Jaccard指数相比，Sørensen相似性指数（Sørensen similarity index）认为两个样方之间共享的物种数量更重要，因此它计算两次。计算公式如下：

转换为Sørensen相异度（Sørensen dissimilarity）：

Simpson相似性/相异度

在两个样方的物种丰富度指数差异很大的情况下（即一个样方比另一个样方具有更多的物种），Simpson相似性指数（Simpson similarity index）更为适用。在这种情形下，如果使用Jaccard或Sørensen相似性指数，它们的值通常非常低，因为会出现分母过大的情况（具有很多的非共享物种，特别是高物种丰富度的样方所贡献）导致指数的总值过低。Simpson相似性指数通过从分数b和c中仅取较小的数据来消除这个问题。计算公式如下：

转换为Simpson相异度（Simpson dissimilarity）：

注意：这里的Simpson相似性指数（或Simpson相异度），不同于Alpha多样性指数中的Simpson指数。

定量非对称相似性指数

例如相似百分率（Percentage similarity），由“1 - 相异百分率”获得（即直接通过“S = 1-D”等转化），相异百分率又称Bray-curtis距离，详见下文“Bray-curtis距离”。

定量对称距离测度

欧几里得距离

欧几里得距离（Euclidean distance）属于欧式距离的一种，是多变量分析中经常使用的一种距离测度，如在线性排序方法PCA、RDA，以及某些层次聚类算法中。欧几里得距离没有上限，最大值取决于数据。

欧几里得距离计算公式如下：

其中，y_1j和y_2j分别是对象1和2中元素j的数值，p是总元素个数。若对应了物种数据，那么y_1j和y_2j即分别是样方1和2中物种j的丰度，p是物种数（样方-物种矩阵中的物种数）。如下展示了仅包含两个物种的两个群落之间的欧几里得距离的计算过程。

然而事实上，在物种数据的分析中，欧几里得距离通常表现得不很理想。主要原因在于它是一个对称的指数，不能有效地处理“双零”问题，而是将“双零”现象视作相同存在的方式处理。当两个样方之间存在很多同为0丰度的物种的情况时，将会很糟糕。关于“双零”以及“对称指数和非对称指数”，可参考前文。

另外一方面，欧几里得距离还具有对“物种丰度的差异”比对“物种是否存在”更加敏感的属性。因此，在使用欧几里得距离评估群落相似性时，两个共享很少物种的样方可能看起来比两个共享较多物种但物种的丰度差别很大的样方具有更相似（具有更低的欧几里得距离），干扰我们对群落相似度的判断。

我们举一个例子说明吧。在下面的物种组成矩阵中，样方1和2不共享任何物种；样方1和3共享物种2、3但丰度不同，物种1均为0。然后我们根据欧几里得距离的公式计算群落距离。

根据欧几里得距离，我们得到了这样的结果：样方1和2具有比样方1和3更高的群落相似度。这显然是很难让人接受的，毕竟样方1和2没有共享任何的物种，而样方1和3共享所有物种，仅仅是丰度相差较大而已。

但如果仍要将欧几里得距离应用在物种数据的分析中，常见的解决方法是首先对原始物种数据执行预转化（常见的如弦转化、Hellinger转化等），然后再使用转化后的数据计算欧几里得距离（即对应于下文提到的弦距离、Hellinger距离，它们仍然是欧氏距离的属性）。尽管弦距离、Hellinger距离等然是对称指数的范畴，但是相较于使用原始物种丰度数据所得的欧几里得距离，弦距离、Hellinger距离的优势体现在存在距离的“上限”，降低了欧几里得距离对“物种丰度”的敏感性，有效减少了“双零”问题导致的误差。

更多情况下，我们在处理物种数据时，会尽可能避开使用欧几里得距离这类的对称指数。例如，通常我们选择使用非对称的Bray-curtis距离等，详见下文。除非特定需要，不得不使用欧氏距离的情况下，可再考虑先转化数据再求欧几里得距离的方法。

但欧几里得距离仍然广泛应用于生态学数据分析中。毕竟对于非物种数据类型，欧几里得距离则通常是广泛适用的，例如，地理因素、环境指标、生物性状等。它们不同于物种，“0值”确实代表了“不存在”这种状态。因此，欧几里得距离可能是多变量分析中最受欢迎的距离测度之一。

此外，对于欧几里得距离（以及下文即将提到的弦距离、Hellinger距离、卡方距离），由于这些距离类型的最大取值远大于1，无法直接通过“S = 1-D”等形式转换为相似度。

因此，可以首先对所有距离指数统一执行标准化处理，例如：

D_norm = D/D_max，或D_norm = (D-D_min)/(D_max-D_min)

此时标准化后的距离矩阵中，所有取值范围将落在[0,1]范围内，就可以转化为相似性了：

S = 1- D_norm，或S = (1- D_norm)^1/2，或S = (1- D²_norm )^1/2

同样地，详见前文的阐述。

弦距离

弦距离（Chord distance）也是欧式距离的一种类型，是根据范数标准化的物种数据计算的欧几里得距离（首先对原始物种丰度数据执行范数标准化，又称为弦转化，再使用弦转化后的物种数据计算欧几里得距离，即为弦距离）。

弦转化意味着多维空间中的物种向量具有单位长度（对于每个样方，物种丰度平方和为1）。相较于使用原始物种数据直接计算的欧几里得距离（没有上限），弦距离具有上限（上限2^0.5）。

弦距离公式如下：

其中，y_1j和y_2j即分别是样方1和2中物种j的丰度，p是物种数（样方-物种矩阵中的物种数）。

Hellinger距离

Hellinger距离（Hellinger distance）也是欧式距离的一种类型，是指通过Hellinger转化后物种数据计算的欧几里得距离（首先对原始物种丰度数据执行Hellinger转化，再使用转化后的物种数据计算欧几里得距离）。相较于使用原始物种数据直接计算的欧几里得距离（没有上限），Hellinger距离具有上限（上限2^1/2）。

对于样方i中物种j，执行Hellinger转化公式如下：

其中，y_ij是样方i中物种j的丰度，y_i+是样方i中所有物种的丰度之和，y^’_ij是样方i中物种j的Hellinger转化后的丰度。从公式可以清楚地看出，它消除了样方间绝对丰度的差异（标准化为相对丰度），且平方根降低了优势物种的重要性。

之后，再基于Hellinger转化后的数据计算的欧几里得距离，即为Hellinger距离。

因此Hellinger距离公式可直接写为：

其中，y_1j和y_2j分别是样方1和2中物种j的丰度，p是物种数（样方-物种矩阵中的物种数）。

Hellinger转化的另一种说法是平方根转化后的弦转化，反过来说，弦转化是多度数据平方后的Hellinger转化。这种关系表明弦距离和Hellinger距离存在关联。

Hellinger转化是一种预先转化物种组成数据以用于线性排序方法的方法，并且被认为是处理带有很多零值的生态数据的合适方法之一。例如，Hellinger距离常应用于物种数据的PCA、RDA等分析中（tb-PCA、tb-RDA），可有效避免“马蹄形效应”的产生。

卡方距离

卡方距离（chi-square distance）通常应用于单峰模型的排序方法中，如CA、CCA。卡方距离没有上限，最大值取决于数据。

卡方距离计算计算公式：

其中，y_1j和y_2j即分别是样方1和2中物种j的丰度；y_+j是物种j在所有样方中的总丰度；y₁₊是样方1中所有物种的总丰度，y₂₊是样方2中所有物种的总丰度，y₊₊是所有样方中所有物种的总丰度；p是物种数（样方-物种矩阵中的物种数）。

定量非对称距离测度

Bray-curtis距离

Bray-curtis距离（Bray-curtis distance）或称Bray-curtis相异度（Bray-curtis dissimilarity）、相异百分率（Percentage difference）。在群落分析中，对于物种数据而言，Bray-curtis距离可能是使用的最广的一种距离测度了。

其计算公式如下：

其中p是物种数（样方-物种矩阵中的物种数），y_ij和y_ik表示两个样方中对应的物种多度。

Bray-curtis距离的取值范围范围由0（两个群落的物种类型和丰度完全一致）到1（两个群落不共享任何物种），因此它也可以直接通过“S = 1-D”等转化为相似性指数（上文提到的“相似百分率”）。Bray-curtis距离适用于群落物种数据分析的原因在于它是一个非对称指数，可有效忽略双零。

这里我们可以再使用Bray-curtis距离重复上文中欧几里得距离中的那个例子。过程就省略了，直接看所得距离测度结果：样方1vs样方2：1，样方2vs样方3：1，样方1vs样方3：0.714，即不共享任何物种的样方1和2、样方2和3的距离达到最大值，而存在共享物种但丰度不同的样方1和3的距离就小于前两者了。相较于上文欧几里得示例，这里的结论显然是更普遍被接受的。

Unifrac距离

这里再额外简介一种特殊的距离，Unifrac距离（Unifrac distance）（Catherine and Rob，2005），它常用于微生物群落的研究中（例如，16S扩增子测序）。上述距离测度的计算方法，仅考虑了物种的存在与否及其丰度，没有考虑物种之间的进化关系，距离0表示两个群落的物种组成结构完全一致。在Unifrac距离中，除了关注考虑了物种的存在与否及其丰度外，还将物种之间的进化关系考虑在内，距离0更侧重于表示两个群落的进化分类完全一致。

例如在16S扩增子测序中，根据16S序列组成构建OTUs进化树，OTUs之间存在进化上的联系，因此不同OTUs之间的（系统发育）距离实际上有“远近”之分。将系统发育树和OTUs丰度数据一起添加至距离的计算中，计算Unifrac距离。而若使用上述提到的其它只基于OTUs丰度数据计算群落距离的方法，则忽略了OTUs之间的进化关系，认为OTUs间的关系平等。当然，并不是说Unifrac距离是最合适16S群落分析的，很多情况下它可能也并没有比只基于OTUs丰度数据计算群落距离的方法（如Bray-curtis距离）“更好”，总之具体问题具体分析吧，根据实际情况选择合适的距离测度。

Unifrac距离分为非加权Unifrac距离（Unweighted unifrac distance）和加权Unifrac距离（Weighted unifrac distance）。两种的主要区别是否考虑了物种的丰度。非加权Unifrac距离只考虑了物种有无的变化，不关注物种丰度，若两个微生物群落间存在的物种种类完全一致，则距离为0；加权Unifrac距离同时考虑物种有无和物种丰度的变化，若两个微生物群落间存在的物种种类及丰度完全一致，则距离为0。

对于非加权和加权Unifrac距离的选择，看网上很多帖子给的经验性建议：在环境样本的检测中，由于影响因素复杂，群落间物种的组成差异更为剧烈，因此往往采用非加权方法进行分析。例如，微生物群落地理学的大尺度研究中，在中国不同气候区采样，比较土壤菌群差异。但如果要研究对照与实验处理组之间的关系，例如研究短期青霉素处理后，人肠道的菌落变化情况，由于处理后群落的组成一般不会发生大改变，但群落的丰度可能会发生大变化，因此更适合用加权方法去计算。

当然这些只是建议，实际情况中可能效果并没有那么好，个人体会如此，比方说你使用加权方法比较中国不同气候区土壤菌群的差异或许也能得到比非加权方法方法更明显的结论。总之具体问题具体分析吧。

相似性指数或距离指数的计算

示例文件见百度盘，提取码namh。

https://pan.baidu.com/s/1DqIYWQ7Arq8TQqWS6yOakw

其中“otu_table.txt”为某16S扩增子测序所得OTU丰度表格（近似为物种丰度表），“otu_tree.tre”为使用各OTU代表序列构建的进化树（即OTU水平的16S进化树）。下面使用其中1个或2个文件展示群落相似性或距离在一些常用软件中的计算方法。

R语言计算

一些R包，如stats包、ade4包、vegan包、cluster包、FD包等，提供了一系列计算相似性指数或距离指数的命令。

这里以vegan包为例作简要展示。常见的相似性指数或距离指数，如Jaccard相似性指数、欧几里得距离、Bray-curtis距离等，可通过vegdist()函数计算。但是需要注意一点，像Jaccard指数这些，本身属性是相似指数，但vegdist()函数的输出结果统统为距离指数，必要时需要通过“S = 1-D，或S = (1-D)^1/2，或S = (1-D²)^1/2”等转换；本身属性是相异指数的，则无需再作转换。

##vegan 包 vegdist() 计算群落相似性/距离，详情 ?vegdist
#以下相似度（S）和距离测度（D）的转换公式，均使用 S = 1-D
library(vegan)

#读入物种数据
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t')
otu <- t(otu)

#jaccard 相异度
jacc_dis <- vegdist(otu, method = 'jaccard', binary = TRUE)
#jaccard 相似度
jacc_sim <- 1 - jacc_dis

#欧几里得距离
eucl_dis <- vegdist(otu, method = 'euclidean')
#这种距离转换为相似性时，需要首先作个标准化，例如
eucl_dis_norm <- eucl_dis/max(eucl_dis)
eucl_sim <- 1 - eucl_dis_norm

#弦距离
otu_chord <- decostand(otu, method = 'normalize')
chord_dis <- vegdist(otu_chord, method = 'euclidean')

#Hellinger 距离
otu_hell <- decostand(otu, method = 'hellinger')
hell_dis <- vegdist(otu_hell, method = 'euclidean')

#卡方距离
#先做卡方标准化，再计算欧几里得距离，即可得
otu_chi <- decostand(otu, method = 'chi.square')
otu_chi_dis <- vegdist(otu_chi, 'euclidean')

#Bray-curtis 距离
bray_dis <- vegdist(otu, method = 'bray')
#Bray-curtis 相似度
bray_sim <- 1 - bray_dis

#更多示例不再展示了
#以上结果均为 dis 类型，若要输出在本地，需要首先转化为 matrix 类型
#以 Bray-curtis 距离为例输出距离测度矩阵
bray <- as.matrix(bray_dis)
write.table(bray, 'bray-curtis_distance.txt', col.names = NA, sep = '\t', quote = FALSE)

#距离矩阵热图
heatmap(bray, scale = 'none', RowSideColors = rep(c('red', 'blue'), each = 18), ColSideColors = rep(c('red', 'blue'), each = 18))

本示例的样本重复性其实不是很好，算了就这样吧。

对于特殊的Unifrac距离，在R中可通过phyloseq包UniFrac()函数计算，GUniFrac包GUniFrac()函数也可以。

##phyloseq 包 distance() 计算群落相似性/距离，详情查看 ?distance
library(phyloseq)

#读入物种数据
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t')
#读入进化树文件（这里为有根树）
tree <- read_tree('otu_tree.tre')
#不过师兄说无根树也可以，只是每次计算的 UUF 和 WUF 会不一样，他会随机指定根

#首先构建 phyloseq 对象
physeq <- phyloseq(otu_table(otu, taxa_are_rows = TRUE), phy_tree(tree))

#计算加权 Unifrac 距离
wei_unif_dis <- distance(physeq, method = 'wunifrac')
#计算非加权 Unifrac 距离
unwei_unif_dis <- distance(physeq, method = 'unifrac')

#也可用其计算其它类型距离，如 Bray-curtis 距离
bray_dis <- distance(physeq, method = 'bray')

#更多示例不再展示了
#注：vegan 包 vegdist() 中的 method 参数同样适用 phyloseq 包 distance() 中的 method 参数

#以上结果均为 dis 类型，若要输出在本地，需要首先转化为 matrix 或类型
#以加权 Unifrac 距离为例
wei_unif <- as.matrix(wei_unif_dis)
write.table(wei_unif, 'weighted-unifrac_distance.txt', col.names = NA, sep = '\t', quote = FALSE)

##GUniFrac 包 GUniFrac() 计算微生物群落的 Unifrac 距离，详情查看 ?GUniFrac
library(GUniFrac)

#读入物种数据
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t')
otu <- data.frame(t(otu))
#读入进化树文件（这个只能使用有根树）
tree <- read.tree('otu_tree.tre')

#计算加权 Unifrac 距离
unifrac <- GUniFrac(otu, tree)$unifracs
wei_unif_dis <- as.dist(unifrac[, , 'd_1'])
#计算非加权 Unifrac 距离
unwei_unif_dis <- as.dist(unifrac[, , 'd_UW'])

#以上结果均为 dis 类型，若要输出在本地，需要首先转化为 matrix 或类型
#以非加权 Unifrac 距离为例
unwei_unif <- as.matrix(unwei_unif_dis)
write.table(unwei_unif, 'unweighted-unifrac_distance.txt', col.names = NA, sep = '\t', quote = FALSE)

qiime程序计算

在基于16S扩增子测序等的微生物群落的分析中，常使用到qiime，qiime中的命令“beta_diversity.py”可以计算常用的相异指数。

#qiime 计算群落距离，shell 命令行操作

cp otu_table.txt otu_table.tsv

#在 otu_table.tsv 开头添加一行“# Constructed from biom file”，以便将 otu_table.tsv 转为 qiime 方便识别样式
sed -i '1i\# Constructed from biom file' otu_table.tsv

#otu_table.tsv 转换为 otu_table.biom
biom convert -i otu_table.tsv -o otu_table.biom --table-type="OTU table" --to-json

#beta_diversity.py 计算常用的 weighted unifrac、unweighted unifrac 以及 Bray-Curtis 距离
beta_diversity.py -i otu_table.biom -o result/ -t otu_tree.tre -m bray_curtis,weighted_unifrac,unweighted_unifrac

#得到的结果输出在“result”路径中
#weighted_unifrac_otu_table.txt、unweighted_unifrac_otu_table.txt、bray_curtis_otu_table.txt

相似性或距离的计算中的注意事项

见前文Beta多样性与生态相似性的结尾注意事项部分。

参考资料

张金屯. 数量生态学. 科学出版社, 2004.

DanielBorcard, FranoisGillet, PierreLegendre, et al. 数量生态学:R语言的应用（赖江山 译）. 高等教育出版社, 2014.

David Zelený博士：https://www.davidzeleny.net/anadat-r/doku.php/en:similarity

Jari Oksanen1. Multivariate Analysis in Ecology - Lecture Notes -. 2004

Legendre P, Legendre L. Numerical Ecology. Second English edition. Developments in Environmental Modelling, 1998, 20, Elsevier

Catherine L, Rob K. UniFrac: a New Phylogenetic Method for Comparing Microbial Communities. Appl Environ Microbiol, 2005, 71(12): 8228-8235.

友情链接